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Abstract 

The paper presents a survey of mathematical problems, techniques, 
and challenges arising in the Thermoacoustic (also called Photoacous- 
tic or Optoacoustic) Tomography. 

1 Introduction 

Computerized tomography has had a huge impact on medical diagnostics. 
Numerous methods of tomographic medical imaging have been developed 
and are being developed (e.g., the "standard" X-ray, single-photon emission, 
positron emission, ultrasound, magnetic resonance, electrical impedance, op- 
tical) [55j EH EE EI]- The designers of these modalities strive to in- 
crease the image resolution and contrast, and at the same time to reduce 
the costs and negative health effects of these techniques. However, these 
goals are usually rather contradictory. For instance, some cheap and safe 
methods with good contrast (like optical or electrical impedance tomogra- 
phy) suffer from low resolution, while some high resolution methods (such 
as ultrasound imaging) often do not provide good contrast. Recently re- 
searchers have been developing novel hybrid methods that combine different 
physical types of signals, in hope to alleviate the deficiencies of each of the 
types, while taking advantage of their strengths. The most successful exam- 
ple of such a combination is the Thermoacoustic Tomography (TAT) 
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|±|62J. Albeit not being a common feature in clinics yet, TAT scanners are 
actively researched, developed and already manufactured, for instance by 
OptoSonics, Inc. ( |http://www.optosonics.com/[ ), founded by the pioneer of 
TAT R. Kruger. 

After a substantial effort, major breakthroughs have been achieved in the 
last couple of years in the mathematical modeling of TAT. The aim of this 
article is to survey this recent progress and to describe the relevant models, 
mathematical problems, and reconstruction procedures arising in TAT, and 
to provide references to numerous research publications on this topic. 

The main thrust of this text is toward mathematical methods; consid- 
erations of the text length, as well as authors' background do not let us 
discuss in any detail industrial and physical set-ups and parameters of the 
TAT technique, and limitations of the corresponding mathematical models. 
Fortunately, the excellent recent surveys by M. Xu and L.-H. V. Wang |117] 
and by A. A. Oraevsky and A. A. Karabutov [87, 88J accomplish all of these 
tasks, and thus the reader is advised to consult with them for all such is- 
sues (see also the recent textbook [1 13] ) . On the other hand, in spite of 
the significant recent progress in mathematics of TAT, there is no compre- 
hensive survey text addressing in details the relevant mathematical issues, 
although the surveys [88l I117j do mention some mathematical reconstruction 
techniques. 

The structure of the paper is a follows: Section [2] contains a brief descrip- 
tion of the TAT procedure. The next Section [3] provides the mathematical 
formulation of the TAT problem. In general, it is formulated as an inverse 
problem for the wave equation. However, in the case of the constant sound 
speed, it can be also described in terms of a spherical mean operator (a 
spherical analog of the Radon transform). The section also contains the list 
of natural questions to be addressed concerning this model. These issues are 
addressed then one by one in the following sections. In particular, Section 
H] discusses uniqueness of reconstruction, i.e. the question of whether the 
data collected in TAT is sufficient for recovery of the information of interest. 
Albeit, for all practical purposes this issue is resolved in Corollary [2J we pro- 
vide an additional discussion of unresolved uniqueness problems, which are 
probably of more academic interest. Section [5] addresses inversion formulas 
and algorithms. In Section [6] effects of having only partial data are discussed. 

1 TAT is also called Photoacoustic (PAT) or Optoacoustic (OAT) Tomography and is 
sometimes abbreviated as TCT, which stands for Thermoacoustic Computed Tomography 
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Section [7] contains results concerning the so called range conditions, i.e. the 
conditions that all ideal data must satisfy. Section [S] provides additional 
remarks and discussions of the issues raised in the previous sections. The 
paper ends with an Acknowledgments section and bibliography. Concerning 
the latter, we need to mention that the engineering and biomedical literature 
on TAT is rather vast and no attempt has been made in this text to create 
a comprehensive bibliography of the topic from the engineering prospective. 
The references in [83 EU [T^J [TT21 HZ] to a large extent fill this gap. The 
authors, however, have tried to present a sufficiently complete review of the 
existing literature on mathematics of TAT. 

2 Thermoacoustic tomography 

In TAT, a short duration EM pulse is sent through a biological object (e.g., 
woman's breast in mammography) with the aim of triggering a thermoa- 
coustic response in the tissue. As it is explained in [If 7j . the radiofrequency 
(RF) and the visible light frequency ranges are currently considered to be the 
most suitable for this purpose. Since mathematics works exactly the same 
way in both of these frequency ranges, we will not make such distinction and 
will be talking about just "an EM pulse". E.g., in Figure [1] a microwave 
pulse is assumed. In most cases the pulse is spatially wide, so that the 
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Figure V. The TAT procedure. 

whole object is more or less uniformly irradiated. Some part of EM energy 
is absorbed throughout the object. The amount of energy absorbed at a 
location strongly depends on local biological properties of the cells. Oxygen 
saturation, concentration of hemoglobin, density of the microvascular net- 
work (angiogenesis), ionic conductivity, and water content are among the 
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parameters that influence the absorption strongly |117j . Thus, if the energy 
absorption distribution function f(x) were known, it would provide a great 
diagnostic tool. For instance, it could be useful for detecting cancerous cells 
that absorb several times more energy in the RF range than the healthy 
ones [621 EH H15j I117j . However, as an imaging tool neither RF waves, nor 
visual light alone would provide acceptable resolution. In the RF case, this 
is due to the long wave length. One can use shorter microwaves, but this 
will be at the expense of the penetration depth. In the optical region, the 
problem is with the multiple scattering of light. So, a different mechanism, 
the so called Photoacoustic Effect [HH I107[ 11131 1117] , is used to image f(x). 
Namely, the EM energy absorption results in thermoelastic expansion and 
thus in a pressure wave p(x, t) (an ultrasound signal) that can be measured 
by transducers placed around the object. Now one can attempt to recover the 
function f(x) (the image) from the measured data p(x, t). Such a measuring 
scheme, utilizing two types of waves, brings about the high resolution of the 
ultrasound diagnostics and the high contrast of EM waves. It overcomes the 
adverse effect of the low contrast of ultrasound with respect to soft tissue. 
In fact, such a low contrast is a good thing here, allowing one to assume in 
the first approximation that the sound speed is constant. This often used 
approximation is not always appropriate, but it is the most studied case at 
the moment. Later on in this text we will describe some initial considerations 
of the variable sound speed case, following [U I5T] . 

For this TAT method (and in particular, for the mathematical model 
described below) to work, several conditions must be met. For instance, the 
time duration of the EM pulse must be shorter than the time it takes the 
sound wave to traverse the smallest feature that needs to be reconstructed. 
The ultrasound detector must be able to resolve the time scale of the duration 
of the EM pulse. On the other hand, the transducer must be also able to 
detect much lower frequencies. Thus, one needs to have extra-wide-band 
transducers, and these are currently available. One can find the technical 
discussion of all these issues, for instance, in [HU 1117] . In this text we will 
assume that all these conditions are met and thus the mathematical models 
described are applicable. 

In the next section we present a mathematical description of the relation 
between f(x) and p(x, t) (similar mathematical problems arise in sonar [73] 
and radar [81] imaging, as well as in geophysics [27]). 
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3 Mathematical model of TAT: 

wave equation and the spherical mean trans- 
form 

3.1 The wave equation model 

We assume that the ultrasound speed at location x is equal to c(x). Then, 
modulo some constant coefficients that we will assume all to be equal to 1, 
the pressure wave p(x, t) satisfies the following problem for the standard wave 
equation [281 flUTl fTT5] : 

(p tt = c 2 (x)A xP ,t>0,xeR 3 

<p(x,0) = f(x), (1) 
(pt{x,0) = 

The goal is to find, using the data measured by transducers, the initial value 
f(x) at t — of the solution p(x, t). 

In order to formalize what data is in fact measured, one needs to specify 
what kind of transducers is used, as well as the geometry of the measurement. 
By the geometry of the measurement we mean the distribution of locations 
of transducers used to collect the data. 

We briefly describe here the commonly considered measurement proce- 
dure, which uses point detectors. Line and planar detectors have also been 
suggested (see Section 18. 1 . lj) . It is too early to judge which one of them 
will become most successful, but the one using point transducers has been 
more thoroughly studied mathematically and experimentally, and thus will 
be mostly addressed in this article. In this case, the transducers are assumed 
to be point-like, i.e. of sufficiently small dimension. A transducer at time 
t measures the average pressure over its surface at this time, which for the 
small size of the transducer can be assumed to be just the value of p(y, t) at 
the location y of the transducer. Dimension count shows immediately that in 
order to have enough data for reconstruction of the function f(x), one needs 
to collect data from the transducers' locations y running over a surface S in 
IR 3 . Thus, the data at the experimentalist's disposal is the function g(y,t) 
that coincides with the restriction of p(x,t) to the set of points y G S. 

Taking into account that the measurements produce the values g(y,t) of 
the pressure p(x, t) of |T| on 5 x IR + , the set of equations (JTJ extends to 
become 
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' Ptt = c 2 {x)A x p,t >0,xe M 3 
< p(x,0) = f(x), 
* Pt(x,0) = 

My,t) = g(y,t),yeSxR+ 



(2) 





X 



Figure 2: An illustration to ([2j). 



The problem now reduces to finding the initial value f(x) in (121) from the 
knowledge of the lateral data g(x,t) (see Figure EH]). A person familiar with 
PDEs might suspect first that there is something wrong with this problem, 
since we seem to have insufficient data for the recovery of the solution of the 
wave equation in a cylinder from the lateral values alone. This, however, 
is an illusion, since in fact there is a significant additional restriction: the 
solution holds in the whole space, not just inside the cylinder S x R + . We 
will see soon that in most cases, the data is sufficient for recovery of f(x). 

3.2 Spherical mean model 

We now introduce an alternative formulation of the problem that works in 
the constant speed case only. We will assume that the units are chosen in 
such a way that c(x) = 1. The known Poisson-Kirchhoff formula [251 Ch. 
VI, Section 13.2, Formula (15)] for the solution of ([T]) gives 



p(x,t) 



c- (t(Rf)(x,t)) , 



(3) 



where 
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is the spherical mean operator applied to the function f(x), and dA is the 
normalized area element on the unit sphere in M. 3 . Hence, knowledge of the 
function g(x,t) for x G S and all t > essentially means knowledge of the 
spherical mean Rf(x,t) at all points (x,t) G S x R + . One thus is lead to 
studying the spherical mean operator R : f —> Rf and in particular its 
restriction R$ to the points x G S only (these are the points where we place 
transducers) : 

Rsf(x,t)= J f(x + ty)dA(y),xeS,t>0. (5) 
\y\=i 

This explains why, in many works on TAT, the spherical mean operator has 
been the model of choice. Albeit the (unrestricted) spherical mean operator 
has been studied rather intensively and for a long time (e.g., [TTJ, [251 EH]), 
its version R s with the centers restricted to a subset S appears to have been 
studied since early 1990s only H-HH, [ISlESlEniEIlIMlESlEaESlESlESl 

and offers quite a few new and often hard questions. 

In what follows, we will alternate between these two (PDE and integral 
geometry) interpretations of the TAT model, since each of them has its own 
advantages. 

3.3 Main mathematical problems of TAT 

We now formulate the typical list of problems one would like to address in 
order to implement the TAT reconstruction. 

1. For which sets S G M 3 is the data collected by transducers placed along 
S sufficient for unique reconstruction of /? In terms of the spheri- 
cal mean operator, the question is whether Rs has zero kernel on an 
appropriate class of functions, say continuous with compact supports. 

2. If the data collected from S is sufficient, what are inversion formulas 
and algorithms? 

3. How stable is the inversion? 

4. What happens if the data is "incomplete"? 
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5. What is the space of all possible "ideal" data g(t, y) collected on a sur- 
face SI Mathematically (and in the constant sound speed case) it is 
the question of describing the range of the operator R$ in appropriate 
function spaces. This question might seem to be unusual (for instance, 
to people used to partial differential equations), but in tomography im- 
portance of knowing the range of Radon type transforms is well known. 
Such information is used to improve inversion algorithms, complete in- 
complete data, discover and compensate for certain data errors, etc. 

(e.g., paEsiEsiiiaEaEiiEaizsiizziisa). 

4 Uniqueness of reconstruction 

Many of the problems of interest to TAT can be formulated in any dimension 
d, albeit the practical dimensions are only d = 3 and d = 2. We will consider 
an arbitrary dimension d whenever we see this suitable. 

Let S C M. d be the set of locations of transducers and / be a compactly 
supported function (one can show that for purposes of uniqueness of recon- 
struction problem, one can always assume that / is smooth |7J). Does the 
absence of the signal on the transducers, i.e. g(t,y) = for all t and y in 
S, imply that / = 0? If the answer is a "yes," we call S - a uniqueness 
set, otherwise a non-uniqueness set. In other words, in terms of TAT, the 
uniqueness sets are those that distributing transducers along them provides 
enough data for unique reconstruction of the function f(x). 

In terms of the wave equation, uniqueness sets are the sets of complete 
observability, i.e. such that observing the motion on this set only, one gets 
enough information to reconstruct the whole oscillation. In terms of the 
spherical mean operator, the question is of whether the equality Rgf = 
implies that / = 0. 

We will address this problem for the constant sound speed case first. 

4.1 Constant speed case 

As it has been discussed, the dimension count makes it clear that S must 
be (d — l)-dimensional, i.e. a surface in 3D or a curve in 2D. We will 
also see that most of such surfaces are "good", i.e. are uniqueness ones (or, 
in other words, provide enough information for reconstruction). Thus, we 
should rather discuss the problem of describing the "bad", non-uniqueness 
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sets. The following simple statement is very important and not immediately 
obvious. 



Lemma 1 |7j [77], [7|] \121f Any non-uniqueness set S is a set of zeros of a 



(non-trivial) harmonic polynomial. In particular, 

1. If there is no non-zero polynomial vanishing on S, then S is a unique- 
ness set. 

2. If there is no non-zero harmonic function vanishing on S, then S is a 
uniqueness set. 

The proof of this lemma is very simple. It works under the assumption 
of exponential decay of the function f(x), not necessarily of compactness of 
its support. It also introduces some polynomials that play significant role in 
the whole analysis of the spherical mean operator R$. 

Let k > be an integer. Consider the convolution 

Q k (x) = \x\ 2k * f(x) = J \x - y\ 2k f(y)dy. (6) 

This is clearly a polynomial of degree at most 2k. Rewriting the integral in 
polar coordinates centered at x and using radiality of \x — y\, one sees that 
Qk{x) is determined if we know the values Rf(x,t) of the spherical mean of 
/ centered at x: 



Q k (x) = c d J t^+^Rfix^dt. 



In particular, If R$f = 0, then each polynomial Qk vanishes on 5*. 

Another observation that is easy to justify is that if the function / is 
exponentially decaying (e.g., is compactly supported), then if all polynomials 
Qk vanish identically, the function itself must be equal to zero. (This is not 
necessarily true anymore if / and its derivatives decay only faster than any 
power, rather than exponentially.) 

Thus, we conclude that if / is not identically equal to zero, then there is at 
least one non-zero polynomial Qk- Since, as we discussed, equality Rsf = 
implies that Qk\s = 0, we conclude that S must be algebraic. 

Now notice the following simple to verify equality (with a non-zero con- 
stant cj(.): 

AQfc = c k Qk-i, (7) 
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where A is the Laplace operator. This implies that the lowest k non-zero 
polynomial Qk is harmonic. Since Qk\s — 0, this proves the lemma. 

Consider now the case when S is a closed (hyper-) surface (i.e., the bound- 
ary of a bounded domain). Since, as it is well known, there is no non-zero 
harmonic function in the domain that would vanish at the boundary (the 
spectrum of the Dirichlet Laplace operator is strictly positive), we conclude 
that such 5* is a uniqueness set for harmonic polynomials. Thus, we get the 
following important 

Corollary 2 |7| \6c^ Any closed surface is uniqueness set for the spherical 
mean Radon transform. 

An older alternative proof of this corollary provides an additional insight 
into the problem. We thus sketch it here. Let us assume for simplicity 
that the dimension d > 3 is odd (even dimensions require a little bit more 
work). Suppose that the closed surface S remains stationary (nodal) for the 
oscillation described by ([I]). Since the oscillation is unconstrained and the 
initial perturbation is compactly supported, after a finite time, the interior 
of S will become stationary. On the other hand, we can think that S is 
fixed (since it is not moving anyway). Then, the energy inside 5* must stay 
constant. This is the contradiction that proves the statement of Corollary [2j 

We will see in the next Section that the same method works in some cases 
of variable sound speed, providing the needed uniqueness of reconstruction 
result. 

This corollary resolves the uniqueness problems for most practically used 
geometries. It fails, however, if / does not decay sufficiently fast (see [3], 
where it is shown in which L p (M. d ) classes of functions f(x) closed surfaces 
remain uniqueness sets). 

It also provides uniqueness for some "limited data" problems. For in- 
stance, if S is an open (even tiny) piece of an analytic closed surface E, it 
suffices. Indeed, if it did not, then it would be a part of an algebraic non- 
uniqueness surface. Uniqueness of analytic continuation would show then 
that the whole X is a non-uniqueness set, which we know to be incorrect. 
This result, however, does not say that it would be practical to reconstruct 
using observations from a tiny S. We will see later that this would not lead 
to a satisfactory reconstructions, due to instabilities. 

A geometry sometimes used is the planar one, i.e. detectors are placed 
along a plane S (line in the 2D). In this case, there is no uniqueness of 
reconstruction when the sound speed is constant. Indeed, if f(x) is odd with 
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respect to S, then clearly all measured data g(t, y) will vanish. However, it is 
well known [251 EE! that functions even with respect to S can be recovered. 
What saves the day in TAT is that the object to be imaged is located on one 
side of S. Then, extending f(x) as an even function with respect to S, one 
can still recover it from the data. 

Although, for all practical purposes the uniqueness of reconstruction prob- 
lem is essentially resolved by the Corollary [2J the complete understanding of 
uniqueness problem has not been achieved yet. Thus, we include below some 
known theoretical results and open problems. 

4.1.1 Non-uniqueness sets in M. 2 . 

In this Section, we follow the results and exposition of [7J EH E2] in discussing 
uniqueness sets in 2D. What are simple examples of non-uniqueness sets? As 
we have already mentioned, any line S (or a hyperplane in higher dimensions) 
is a non-uniqueness set, since any function / odd with respect to S will clearly 
produce no signal: R$f = 0. Analogously, consider a Coxeter system T, N of 
N lines passing through a point and forming equal angles (see Fig. [3]). 




Figure 3: Coxeter cross Tj N . 

Choosing the intersection point as the pole and expanding functions into 
Fourier series with respect to the polar angle, it is easy to discover existence 
of an infinite dimensional space of functions that are odd with respect to 
each of the N lines. Thus, such a cross Eat is also a non- uniqueness set. Less 
obviously, one can use the infinite dimensional freedom just mentioned to 
add any finite set $ of points still preserving non-uniqueness. The following 
major and very non-trivial result was conjectured in [7T1 [72] and proven in 
[7J. It shows that there are no other bad sets S besides the ones we have just 
discovered: 



11 



Theorem 3 A set S C M 2 is a non-uniqueness set for the spherical mean 
transform in the space of compactly supported functions, if and only if 

S C ujY, n U $, 

where is a Coxeter system of lines, u is a rigid motion of the plane, and 
$ is a finite set. 

A sketch of a rather intricate proof of this result is provided in Section 

IO 

4.1.2 Higher dimensions 

Here we present a believable conjecture of how the result should look like in 
higher dimensions. 

Conjecture 4 FflA set S C M. d is a non-uniqueness set if and only if S C 
co>X U $, where E is the surface of zeros of a homogeneous harmonic polyno- 
mial, uj is a rigid motion of M. d , and $ is an algebraic surface of dimension 
at most d — 2. 




Figure 4: A picture of a 3-dimensional non-uniqueness set. 

The progress towards proving this conjecture has been slow, albeit some 
partial cases have been treated (PQ-[I2])- E.g., in some cases one can prove 
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that S is a ruled surface (i.e., consists of lines), but proving that these lines 
(rules) pass through a common point remains a challenge. It is known, 
though, that both the zero sets of homogeneous harmonic polynomials and 
algebraic subsets of dimension at most d — 2 are non- uniqueness sets [21 [7] , 
and thus one should avoid using them as placements of transducers for TAT. 

4.1.3 Relations to other areas of analysis 

The problem of injectivity of Rs has relations to a wide variety of areas of 
analysis (see [HE] for many examples). In particular, the following interpre- 
tation is important: 

Theorem 5 |7| The following statements are equivalent: 

1. S cM. d is a non-uniqueness set for the spherical mean operator. 

2. S is a nodal set for the wave equation, i.e. there exists a non-zero 
compactly supported f such that the solution of the wave propagation 



vanishes on S for any moment of time. 

3. S is a nodal set for the heat equation, i.e. there exists a non-zero 
compactly supported f such that the solution of the problem 



vanishes on S for any moment of time. 

The interpretation in terms of the wave equation provides important PDE 
tools and insights, which have lead to a recent progress [331 [12] (albeit it 
has not lead yet to a complete alternative proof of Theorem [3]). The rough 
idea, originally introduced in [33], is that if S is a nodal set, then it might be 
considered as the fixed boundary. In this case, the signals must go around S. 



problem 




u(x,0) = 0, 
Ut(x,0) = f(x) 
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However, in fact, there is no obstacle, so signals can propagate along straight 
lines. Thus, in order to avoid discrepancies in arrival times, S must be very 
special. One can find details in [33] and in [T2] . 

4.2 Uniqueness in the case of a variable sound speed 

It is shown in [351 Theorem 4] that uniqueness of reconstruction also holds in 
the case of a smoothly varying (strictly positive) sound speed, if the source 
function f(x) is completely surrounded by the observation surface S (in other 
words, if there is no US signal coming from outside of S). The proof uses 
the celebrated unique continuation result by D. Tataru [108J. 

One can also establish uniqueness of reconstruction in the case of the 
source not necessarily completely surrounded by S. However, here we need 
to impose an additional non-trapping condition on the sound speed. We 
assume that the sound speed is strictly positive c(x) > c > and such that 
c(x) — 1 has compact support, i.e. c(x) = 1 for large x. 

Consider the Hamiltonian system in Rj^ with the Hamiltonian H = 



The solutions of this system are called bicharacteristics and their projections 
into R™ are rays. 

We will assume that the following non-trapping condition holds: 
all rays (with £ ^ 0) tend to infinity when t — » oo. 

Theorem 6 j^j Under the non-trapping conditions formulated above, com- 
pactly supported function f(x) is uniquely determined by the data g measured 
on S for all times. (No assumption of f being supported inside S is imposed.) 

One should mention that ray trapping can occur for some sound speed 
profiles. For instance, if c(x) = \x\ for some range r\ < \x\ < r2, then there 
are rays trapped in this spherical shell. We are not sure what happens in 
this case to the uniqueness of reconstruction statement of Theorem [6] and 
inversion formula of Theorem [7J 



c 2 W |(M2. 



2 ISI ' 
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5 Reconstruction: formulas and examples 



Here we will address the procedures of actual reconstruction of the source 
f(x) from the data g(t,y) measured by transducers. 

5.1 Constant sound speed 

We assume here that the sound speed is constant and normalized to be equal 
to 1. 

5.1.1 Inversion formulas 

Before we move to our case of interest, which is spheres centered on a closed 
surface 5* surrounding the object to be imaged, we briefly refer to related 
but somewhat different work. Namely, the problem of recovering functions 
from integrals over spheres centered on a (hyper)plane S has attracted a 
lot of attention over the years. Albeit, as it has been mentioned before, 
there is no uniqueness in this case (functions odd with respect to 5* are 
annihilated), even functions can be recovered. Thus also functions supported 
on one side of the plane can be as well, by means of their even extension. 
Many explicit inversion formulas and procedures have been obtained for this 



situation [IHl [26l EH [391 SH EQl [771 EQl [89l [90l HUH QUI US]. We will not 



provide any details here, since this acquisition geometry is not very useful. 
In particular, this is due to "invisibility" of some parts of the interfaces, 
see Section [HI which arises from truncating the plane. The same problem 
is encountered with some other unbounded acquisition surfaces, such as a 
surface of an "infinitely" long cylinder. 

Thus, it is more practical to place transducers along a closed surface 
surrounding the object. The simplest surface of this type is a sphere. 

5.1.2 Fourier expansion methods 

Let us assume that S is the unit sphere in M. n . We would like to reconstruct 
a function f(x) supported inside S from the known values of its spherical 
integrals g(y, r) with the centers on S: 




yeS. 
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The first inversion procedures for the case of spherical acquisition were de- 
scribed in [22] in 2D and in [S3] in 3D. These solutions were obtained by 
harmonic decomposition of the measured data and the sought function, and 
by equating coefficients of the corresponding Fourier series. 

In particular, the 2-D algorithm of [82] is based on the Fourier decompo- 
sition of / and g in angular variables: 

oo 

f(x) = Y J fk(p)e lk ^ x = (pcos(^),psim») (9) 



9(y(0),r) = J29m(r)e lke , y = (Rcos(9),Rsm(9)). 

— oo 

Following [82] we consider the Hankel transform g m ,j(X) of the Fourier coef- 
ficients g m (r) (divided by 27rr) 

ftn,j(A) = J™ g m {r)U\r)dr = H (j^ty ■ (10) 

To simplify the presentation we introduce the convolution Gj(X,y) of the 
sought function with the Bessel function Jq(X\x — y\). 

Gj(X,y)= I f(x)J (X\x-y\)dx, (11) 



One can notice that g m j(X) are the Fourier coefficients of Gj(X,y) in 9: 

9m,jW = o" - r G ^ vY~ m6 d9. (12) 
27r Jo 

Now coefficients f m (p) can be recovered from g m (r) by application of the 
addition theorem for the Bessel function Jo(A|x — y\): 

oo 

J (A|x -y\) = J2 Jm{m)J m {X\y\)e- im ^ e \ (13) 



Indeed, by substituting equations (Ej) and ([TBI into lull) , and (fill) into (fl2 
one obtains 

9m,j{X) = 27rJ m (A|i?|) / f m {p)J m {Xp)pdp = H m {f m {p)), 
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where H m is the m-th order Hankel transform. Since the latter transform 
is self-invertible, the coefficients f m (p) can be recovered by the following 
formula 



fm(p) = ti r , 



9m,j{X) 

[JM\R\) 



:7~to 



9m{r) 
2nr 



(14) 



J m {\\R\Y 

which is the main result of [52] • Function f(x) can now be reconstructed by 
summing series (Q. 

Note that the above method requires a division of the Hankel transform 
of the measured data by Bessel functions J m that have infinitely many zeros. 
Theoretically, there is no problem; the Hankel transform 7i has to 

have zeros that would cancel those in the denominator. However, since the 
measured data always contain some error, the exact cancelation is not likely 
to happen, and one needs a sophisticated regularization scheme to keep the 
total error bounded. 

This problem can be avoided by replacing in (fTUI) Bessel function J by 



Hankel function H, 



(i). 



<Wf(A) 



The addition theorem for Hq \\\x 



2R 



{r)H^\\r)dr. 



y\) takes form 

oo 

# (1) (A|* - y\) = J2 J m (M^\)H^(My\)e- im{ip - d \ 



and by proceeding as before one can obtain the following formula for f m (p): 



fm(p) = % 



H&\\\R\) 



TCr 



H$(\\R\) 



2R 



g m (r)H^\Xr)dr 



Unlike J m , Hankel functions Hm\t) do not have zeros for all real values of t 
and therefore problems with division by zeros do not arise in this amended 
version of the method |82j . 

This derivation can be repeated in 3-D, with the exponentials e lke replaced 
by the spherical harmonics, and with cylindrical Bessel functions replaced by 
their spherical counterparts. By doing this one will arrive at the Fourier 
series method of E3I- Our use of Hankel function above is similar to 



the way the authors of [83] utilized spherical Hankel function h 
the divisions by zero. 



(i) 



to avoid 
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5.1.3 Filtered backprojection methods 

The favorite way of inverting Radon transform for tomography purposes is 
by using filtered backprojection type formulas, which involve filtration in 
Fourier domain followed (or preceded) by a backprojection. In the case of 
the set of spheres centered on a closed surface (e.g., sphere) S, one expects 
such a formula to involve a filtration with respect to the radius variable and 
then some integration over the set of spheres passing through the point of 
interest. For quite a while, no such type formula had been discovered. This 
did not prevent practitioners from reconstructions, since good approximate 
inversion formulas (parametrices) could be developed, followed by an iterative 
improvement of the reconstruction, see e.g. reconstruction procedures in 

[tm ma ebi ma nag , and especially una Eg . 

The first set of exact inversion formulas of the filtered backprojection type 
was discovered in [33]. These formulas were obtained only in odd dimensions. 
Several different variations of such formulas (different in terms of exact order 
of the filtration and backprojection steps) were developed. Let us denote by 
g(p,r) = r 2 Rgf the spherical integral, rather than the average, of /. Then 
various versions of the 3D inversion formulas that reconstruct a function f(x) 
supported inside S from its the spherical mean data Rsf, read: 



f( x ) = -8^r A * I 9{y, \V ~ x\)dA(y), 

dB 




Recently, analogous formulas were obtained for even dimensions in [32]. De- 
noting by g, as before the spherical integrals (rather than averages) of /, the 
formulas of [32J in 2D look as follows: 

2R 

f{x) = 2^R A J J 9(y,t)\og(t 2 -\x-y\ 2 )dtdl(y), (16) 

dB 

or 

2R 

= m II I ~ '* - ^ dt (17 > 

dB 
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A different set of explicit inversion formulas that work in arbitrary dimensions 
was presented in [59"] . 



/0) = 4 ( 27r )n-i div J n (y) h (yi \ x ~ y\) dA (y)- 

dB 



(18) 



Here 



h{y,t) 



/2R 

J Y{\t) ( / J{\t')g(y,t')dt' 
Vo 

"2R 

-J(Xt) ( J Y(\t')g(y,t')dt' ] A 2 "- : \/A. 

.0 



(19) 



J(t) 



J, 



n/2-l 



Y(t) 



Y 



n/2-l 



(0 



pi/2-1 ' v I fn/2-1 ' 

^n/2-i(*) an d Y n /2~i{t) are respectively the Bessel and Neumann functions 
of order n/2 — 1, and n(y) is the vector of exterior normal to dB. 

In 2-D equations (fl"8l . (Tl9|) can be simplified to yield the following recon- 
struction formula 



/(*) 



— ^r-div / n(y) 
2tt 2 J KyJ 



2R 



\x — y\ 2 — t' 2 



dt' 



dl(y). 



A similar simplification is also possible in 3D resulting in the formula 



/(*) 



^divjn(y) 



1 d g(y,t) 
t dt t 



dB 



dA{y). 



(20) 



t=\y-x\ 



Equation (!20j) is equivalent to one of the formulas derived in [116J for the 3D 
case. It is interesting to notice that the "universal" formula of |116j holds for 
all geometries when the backprojection type formulas are known: spherical, 
cylindrical, and planar. It is not very likely that such explicit formulas would 
be available for any closed surfaces S different from spheres (see a related 
discussion in pT5l EI]). 
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5.1.4 Series solutions for arbitrary geometries 

Although, as we have just mentioned, we do not expect such explicit formulas 
to be derived for non-spherical closed surfaces S, there is, however, a different 
approach [70] that theoretically works for any closed S and that is practically 
useful in some non-spherical geometries. 

Let and u m (x) be the eigenvalues and normalized eigenf unctions of 
the Dirichlet Laplacian —A on the interior Q of a closed surface S: 

Au m {x) + \ 2 m u m {x) = 0, xett, flCR", (21) 
u m (x) =0, x G S, 

1 1 Urn 1 1 9 = 



-m. 1 1 2 — J \Urn{p^)\ dx — 1. 

h 



As before, we would like to reconstruct a compactly supported function f(x) 
from the known values of its spherical integrals g{y,r) with the centers on S: 

g{y, r)=//(» + Mr-^, yes. 

We notice that u m (x) is the solution of the Dirichlet problem for the Helmholtz 
equation with zero boundary conditions and the wave number A m , and thus 
it admits the Helmholtz representation 

f d 

u m (x)= $x m (\x -y\)jr; u m {y)ds{y) x E fi, (22) 

J 80, ° n 

where Q\ m (\x —y\) is a free-space rotationally invariant Green's function of 
the Helmholtz equation fl2Tl) . 

The eigenfunctions {u m (x)}'^ ) form an orthonormal basis in L 2 (fi). There- 
fore, f(x) can be represented by the series 

oo 

f{x) = S ^a m u m (x) (23) 

m=0 

with 

a m = / u m (x)f(x)dx. 



n 
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Since f(x) is Cq, series (I23p converges pointwise. A reconstruction formula of 
a m , and thus of f(x), will result if we substitute representation (12"2"1) into (1231) 
and interchange the order of integrations. Indeed, after a brief calculation 
we will get 

ot m = I u m (x)f(x)dx = I I(y,X m )-^u m (y)dA(x), (24) 



n Jan 



where 



Ifo, A) = / $ x (\x-y\)f(x)dx. (25) 

Certainly, the need to know the spectrum and eigenfunctions of the 
Dirichlet Laplacian imposes a severe constraint on the surface S. However, 
there are simple cases when the eigenfunctions are well known, and fast sum- 
mation formulas for the corresponding series are available. Such is the case 
of a cubic measuring surface S (see [70J ) ; the eigenfunctions products 
of sine functions 

8 . vrmixi . 7im 2 x 2 . nm 3 x 3 
u m {x) = sm — ^— sm — ^— sm - , (26) 

where m = (m 1 ,m 2 ,m 3 ), m 1 ,m 2 ,m3 G N, and the eigenvalues are easily 
found as well 

A m = TT 2 \m\ 2 /R 2 . (27) 

Sum (12"51) is just a regular 3-D Fourier sine series easily computable by ap- 
plication of the Fast Sine Fourier transform algorithm. The algorithmic 
trick that allows one to calculate fast the coefficients a m consists in com- 
puting first integrals (125]) on a uniform mesh in A. This is easily done by 
a one-dimensional Fast Cosine Fourier transform algorithm, with <$>\(t) = 
cos(Ai)/t. The normal derivatives of u m (x) are also products of sine func- 
tions, this time two-dimensional ones. This, in turn, permits rapid evaluation 
of integrals J m I(y, \)J^u m (y)dA(x) for each mesh value of A, and for each 
one of the six faces dQi, % — 1, 6 of the cube. Finally, the computation of 
a m using equation (1241) reduces to the interpolation in the spectral param- 
eter A, since the integrals in the right hand side of this equation have been 
computed for the mesh values of this parameter (not for A m ). Due to oscil- 
latory nature of the integrals (125]) a low order interpolation here would lead 
to inaccurate reconstructions. Luckily, however, these integrals are analytic 
functions of parameter A (due to the finite support of g). Hence, high order 
polynomial interpolation is applicable, and numerics yields very good results. 
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The algorithm we just described requires 0(m 3 logm) floating point op- 
erations if the reconstruction is to be performed on an m x m x m Cartesian 
grid, from comparably discretized data measured on a cubic surface. In prac- 
tical terms, it yields reconstructions in the matter of several seconds on grids 
with total number of nods exceeding a million [70] . 

5.1.5 Time reversal (backpropagation) methods 

In the constant speed case, the following approach is possible in 3-D: due 
to the validity of the Huygens' principle (i.e., the signal escapes from any 
bounded domain in finite time), the pressure p(t,x) inside S will become 
equal to zero for any time T larger than the time required to cross the domain 
(i.e., time that it takes the sound to move along the diameter of S, which 
for c = 1 equals the diameter). Thus, one can impose the zero conditions 
on p(t, x) for t = T and solve the wave equation §2§ back in time, using the 
measured data g as the boundary values. The solution of this well posed 
problem at t = gives the desired source function f(x). Such methods have 
been successfully implemented [22] . 

Although in 2D or in presence of sound speed variations, Huygens' princi- 
ple does not hold anymore, and thus the signal theoretically will stay forever, 
one can find good approximate solutions using a similar approach [U [TS] , see 
discussions of the variable speed case below. 

5.1.6 Examples of reconstructions and additional remarks about 
the inversion formulas 

• It is well known that different analytic inversion formulas in tomogra- 
phy can behave differently in numerical implementation (e.g., in terms 
of their stability), However, numerical implementation seems to show 
that the analytic (backprojection type) formulas (IT51) - (|2"0|) . in spite of 
some of them being not equivalent, work equally well. See, for example 
the results of an analytic formula reconstruction in 3D shown in Fig. 

El 

• It is worth noting that although formulas f|T5|) - f|T6|) and ffl8|) -f l20|) will 
yield identical results when applied to functions that can be represented 
as the spherical mean Radon transform of a function supported inside 
S, they are in general not equivalent when applied to functions with 
larger supports. Simple examples (e.g., of / being the characteristic 
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Figure 5: A mathematical phantom in 3D (left) and its reconstruction using 
an analytic inversion formula. 

set of a large ball containing S) show that these two types of formulas 
provide different reconstructions. 

• An interesting observation is that backprojection formulas (I15l)- (l20l) do 
not reconstruct the function / correctly inside the surface S, if / has 
support reaching outside S. For instance, applying the reconstruction 
formulas to the function Rs{x\x\<3) leads to an incorrect reconstruction 
of the value of / = X\x\<3 inside S = {\x\ < 1}. (Here by xv we denote 
the characteristic function of the set V, i.e. it takes the value 1 in V 
and zero outside. So, x\x\<3 is the characteristic function of the ball of 
radius 3 centered at the origin.) 

An another example: if one adds to the phantom shown in Fig. [5] 
two balls to the right of the surrounding sphere S, this leads to strong 
artifacts, as seen on Fig. El 

What is the reason for such a distortion? If one does not know in ad- 
vance that / has support inside S, the backprojection formulas shown 
before use insufficient information to recover a function with a larger 
support, and thus uniqueness of reconstruction is lost. Then the for- 
mulas misinterpret the data, wrongly assuming that they came form a 
function supported inside S and thus reconstructing the function in- 
correctly. 

Notice that the series reconstruction of the preceding Section is free of 
such problem. E.g., the reconstruction shown in Fig. [7| confirms this. 



23 



Figure 6: A perturbed reconstruction, due to presence of two additional balls 
outside S (not shown on the picture). 




Figure 7: In the phantom shown on the left, most disks are located outside 
the square acquisition surface S indicated by the dotted line. This does not 
perturb the reconstruction inside S (right). 

5.2 Reconstruction in the variable speed case 

We will assume here that the sound speed c(x) is smooth, positive, constant 
for large x, and non-trapping. Although most analytic techniques we de- 
scribed above do not work in the variable speed case, some formulas can be 
derived and algorithms can be designed. This work is in a beginning stage 
and the results described below most surely can and will be improved. 

5.2.1 "Analytic" inversions 

Let us denote by Q the interior of the observation surface S, i.e. the area 
where the object to be imaged is located. Consider in Q the operator A = 
—c 2 (x)A with zero Dirichlet conditions on the boundary S = dQ. This 
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operator is self- adjoint, if considered in the weighted space L 2 (Q; c~ 2 (x)). 

We also denote by E the operator of harmonic extension, which trans- 
forms a function on 5* to a harmonic function on Q which coincides with (f> 
on S. 

The following result provides a formula for reconstructing / from the data 

9- 

Theorem 7 ^ The function f(x) in (TJj) can be reconstructed in Q as fol- 
lows: 

oo 

f(x) = (Eg\ t=0 ) - j A-hm(rA^E(g tt )(x,r)dT. (28) 
o 

The validity of this result hinges upon decay estimates for the solution (so 
called local energy decay [29J, I110[ llllj ). which hold under the non-trapping 
condition. These estimates guarantee a qualified decay of the solution p(t, x) 
inside any bounded region, e.g. in Q, when time t increases. In odd dimen- 
sions decay is exponential, but only polynomial in even dimensions. The 
decay can be used instead of Huygens' principle to solve the wave equation 
backwards, starting at the infinite time. This leads to the formula ( 1281) . 

Due to functions of the operator A being involved, it is not that clear how 
explicit this formula can be made. For instance, it would be interesting to 
see whether one can derive from (128]) a backprojection inversion formula for 
the case of a constant sound speed and S being a sphere (we have already 
seen that such formulas are known). 

5.2.2 Backpropagation 

The exponential decay at large values of time can be used as follows: for a 
sufficiently large T, one can assume that the solution is practically zero at 
t = T. Thus, imposing zero initial conditions at t = T and solving in reverse 
time direction, one arrives at t = to an approximation of f(x) [18J. 

5.2.3 Eigenfunction expansions 

One natural way to try to use the formula (f2"8"j) is to use eigenfunction expan- 
sion of the operator A in Q (assuming that such expansion is known). This 
immediately leads to the following result: 
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Theorem 8 Under the same conditions on the sound speed as before, func- 
tion f(x) can be reconstructed inside Q from the data g in (TJ|) ; as the following 
L 2 (B)- convergent series: 

/(z)=^/ fe V> fe (z), (29) 

k 

where the Fourier coefficients fk can be recovered using one of the following 
formulas: 

oo 

fk = Xfg k (0) - X k 3 J sin {\ k t)g'l{t)dt, 
o 

oo 

fk = Xfg k (0) + Xf J cos (X k t)g' k (t)dt, or (30) 
o 

oo oo 

fk = -Xj: 1 f sin (X k t)g k (t)dt = -X]: 1 f f sin (X k t)g(x, t)^(x)dxdt, 

OS 

and 

9k(t) = J g(x,t)^-(x)dx. 
s 

Here v denotes the external normal to S . 

One notices that this is a generalization to the variable sound speed case 
of the expansion method of [70], discussed in Section [5.1.41 An interesting 
feature is that, unlike in [70], we do not need to know the whole space Green's 
function for A (which is certainly not known). 

It is not clear yet how feasible numerical implementation of this approach 
could be. 

6 Partial data. "Visible" and "invisible" sin- 
gularities 

Uniqueness of reconstruction does not imply practical recoverability, since the 
reconstruction procedure might be severely unstable. This is well known to 
be the case, for instance, in incomplete data situations in X-ray tomography, 
and even for complete data problems in some imaging modalities, such as 
the electrical impedance tomography [6U EH [TBI [77] . 
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In order to describe the results below, we need to explain the notion 
of the wave front set WF(f) of a function f(x). This set carries detailed 
information on singularities of f(x). It consists of pairs (x,£) of a point x 
in space and a wave vector (Fourier domain variable) £ 7^ 0. It is easier to 
say what it means that a point (xo, £0) is not in the wave front set WF(f). 
This means that one can smoothly cut-off / to zero at a small distance from 
Xq in such a way that the Fourier transform of the resulting function 

4>(x)f(x) decays faster than any power of £ in directions that are close to 
the direction of £o- We remind the reader that if this Fourier transform 
decays that way in all directions, then f(x) is smooth near the point xq. So, 
the wave front set contains pairs (xo,£o) such that / is not smooth near x > 
and £ indicates why it is not: the Fourier transform does not decay well in 
this direction. For instance, if f(x) consists of two smooth pieces joined non- 
smoothly across a smooth interface E, then WF(f) contains pairs (x, £) such 
that x is in £ and £ is normal to X at x. One can find simple introduction 
to the notions of micro local analysis, such as the wave front set, for instance 
in [TUB] . 

Analysis done in [HS] for the constant speed case (equivalently, for the 
spherical mean transform Rs), showed which parts of the wave front (and 
thus singularities) of a function / can be recovered from its partial X-ray 
data. An analog of this result also holds for the spherical mean transform 
Rs [Z3] (see also |12Uj for a practical discussion). We formulate it below in 
an imprecise form (see [73] for precise formulation). 

Theorem 9 |75]/ A wavefront set point (x,£) of f is "stably recoverable" 
from Rsf if and only if there is a circle (sphere in higher dimensions) centered 
on S , passing through x, and normal to £ at this point. 

As we have already mentioned, this result does not exactly hold the way it is 
formulated and needs to include some precise conditions (see J73J Theorem 
3]). The statement is, for instance, correct if S is a smooth hypersurface and 
the support of / lies on one side of the tangent plane to S at the center of 
the sphere mentioned in the theorem. 

Talking about jump singularities only (i.e., interfaces between smooth 
regions inside the object to be imaged), this result says that in order for a 
piece of the interface to be stably recoverable (dubbed "visible"), one should 
have for each point of this interface, a sphere centered at S and tangent to 
the interface at this point. Otherwise, the interface will be blurred away 
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(even if there is a uniqueness of reconstruction theorem). The reason is that 
if all spheres of integration are transversal to the interface, the integration 
smoothes off the singularity, and therefore its recovery becomes highly unsta- 
ble (numerically, one has to deal with inversion of a matrix with exponentially 
fast decaying singular values) . The Figure [8] below shows an example of an 
incomplete data reconstruction from spherical mean data. One sees clearly 
the effect of disappearance of the parts of the boundaries that are not touched 
tangentially by circles centered at transducers' locations. 




Figure 8: Effect of incomplete data: the phantom (left) and its incomplete 
data reconstruction. The transducers were located along a 180° circular arc 
(the left half of a large circle surrounding the squares). 



7 Range conditions 

As it has already been mentioned, the space of functions g(t,y) that could 
arise as exact data measured by transducers (i.e., the range of the data), 
is very small (of infinite codimension in the spaces of all functions of t > 
0,y G S). Knowing this space (range) is useful for many theoretical and 
practical purposes (reconstruction algorithms, error corrections, incomplete 
data completion, etc.), and thus has attracted a lot of attention (e.g., [30| 
1381 1391 1301 1531 1511 161 1661 1671 1681 171 1761 1771 1781 1901 [TOO]. 
For instance, for the standard Radon transform 

f(x) -»• g(s,u) = J f(x)dx,\uj\ = 1, 

X-LO = S 

the range conditions on g(s,u) are: 



28 



1. evenness: g(—s,—u>) — g(s,u) 

2. moment conditions: for any integer k > 0, the fcth moment 



CO 



Gk(to) — s g(u,s)ds 



extends from the unit circle of vectors wtoa homogeneous polynomial 
of degree k in to. 

The evenness condition is obviously necessary and is kind of "trivial". It 
seems that the only non-trivial conditions are the moment ones. However, 
here the standard Radon transform misleads us, as it often happens. In fact, 
for more general transforms of Radon type it is often easy (or easier) to find 
analogs of the moment conditions, while analogs of the evenness conditions 
are often elusive (see [Ml EH E7J EE [771 El] devoted to the case of SPECT 
(single photon emission tomography)). The same happens in TAT. 

Let us deal first with the case of a constant sound speed, when one can 
think of the spherical mean transform R$ instead of the wave equation model. 
An analog of the moment conditions was already present implicitly (without 
saying that these were range conditions) in [7_H [721 El an d explicitly formu- 
lated as such in [53] • Indeed, our discussion in Section 4 of the polynomials 
Qk provides the following conditions of the moment type: 

Moment conditions |7| [77], [7|| \95f on data g(p, r) = R s f(p, r) look as 



follows: for any integer k > 0, the moment 

oo 

M k {uo) = [ r 2k+d - 1 g(p,r)dr 



can be extended from S to a (non-homogeneous) polynomial Qk{x) of degree 
at most 2k. 

These conditions, however, are incomplete, and in fact infinitely many 
others, which play the role of an analog of evenness, need to be added. 

Complete range descriptions for R$ when S is a circle in 2D were discov- 
ered in [13j and then in odd dimensions in [31]. They were then extended 
to any dimension and interpreted in several different ways in [6|. These 
conditions happen to be intimately related to PDEs and spectral theory. 
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Figure 9: 



In order to describe these conditions, we need to introduce some notations. 
Let B be the unit ball in S - the unit sphere, and C - the cylinder Bx[0,2] 
(see Fig. E>- 

We introduce the spherical mean operator Rs as before: 

Rsf{x,t)= f f(x + ty)dA(y),xeS. 
J\v\=i 

Several different range descriptions for Rs were provided in [6], out of 
which we only show a few: 

Theorem 10 [ffj The following three statements are equivalent: 

1. The function g G C^°(S x [0,2]) is representable as Rsf for some 
f G Cq°(B). (In other words, g represents an ideal (free of errors) set 
of TAT data.) 

2. (a) The moment conditions are satisfied. 

(b) Let —A 2 be any eigenvalue of the Laplace operator in B with zero 
Dirichlet conditions and ip\ be the corresponding eigenf unction. 
Then the following orthogonality condition is satisfied: 

J g(x,t)d u ij x (x)j n/ ^ 1 (\t)t n - 1 dxdt = 0. (31) 

Sx[0,2] 

Here j p (z) = c v ^p- is the so called spherical Bess el function. 

3. (a) The moment conditions are satisfied. 



30 



(b) Let g{x, A) = J g(x,t)j n /2-i{Xt)t n ~ l dt. Then, for any m e Z ; the 
m th spherical harmonic term g m {x, A) of g{x, A) vanishes at all 
zeros A 7^ of Bessel function J m+n /2-i(A). 

Remark 11 

1. In odd dimensions, moment conditions are not necessary, and thus con- 
ditions 2(b) or 3(b) suffice. (A similar earlier result was established for 
a related transform in \3Jfl.) 

2. The range conditions (2) of the previous Theorem are also necessary 
when S is the boundary of any bounded domain, not necessarily a 



3. An analog of these conditions can be derived for a variable sound speed 
(without non-trapping conditions imposed). 

8 Concluding remarks 

8.1 Variations of the TAT procedure 
8.1.1 Planar and linear transducers 

Assuming that transducers are point-like, is clearly an approximation, and 
in fact a transducer measures the average pressure over its area. It has 
been rightfully claimed that the point approximation for transducers should 
lead to some blurring in the reconstructions. This, as well as intricacies of 
reconstructions from the data obtained by point transducers, triggered recent 
proposals for different types of transducers (see [2Ql[2TJ, |37]-[52], [92], [93]). 
In these papers, it was suggested to use either planar, or line detectors. 

In the first case [37], the detectors are assumed to be large and planar, 
ideally assumed to be approximations of infinite planes that are placed tan- 
gentially to a sphere containing the object. Thus, the data one collects is 
the integrals of the pressure over these planes, for all values of t > 0. If one 
takes the standard 3D Radon transform of the pressure p(x, t) with respect 
to x: 



sphere. 
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where dA is the surface measure and u is a unit vector in R 3 , this is well 
known to reduce the 3D Laplace operator A x to the second derivative d 2 /ds 2 
[301 138| [39| BUI 1531 151] , and thus the 3D wave equation to the string vibra- 
tion problem. The measured data provide the boundary conditions for this 
problem. The initial conditions in ([1]) mean evenness with respect to time, 
and thus the standard d'Alambert formula leads to the immediate realiza- 
tion that the measured data is just the 3D Radon transform of f(x). Thus, 
the reconstruction boils down to the well known inversion formulas for the 
Radon transform. 

Another proposal ([201 [21], 02] -[S2], [921 [93]) is to use line detectors 
that provide line integrals of the pressure p(x, t). Such detectors can be 
implemented optically, using either Fabry- Perot [20J, or Mach-Zehnder [95] 
interferometers. 

Suppose that the object is surrounded by a surface that is rotation in- 
variant with respect to the z-axis. It is suggested to place the line detectors 
perpendicular to the z-axis and tangential to the surface. The same consid- 
eration as above then shows that after the 2D Radon (or X-ray, which in 
2D is the same) transform in each plane orthogonal to z-axis, the 3D wave 
equation converts into the 2D one for the Radon data. The measurements 
provide the boundary data. Thus, the reconstruction boils down to solving 
a 2D problem similar to the one in the case of point detectors, and then 
inverting the 2D Radon transform. 

Due to the recent nature of these two projects, it appears to be too early 
to judge which one will be superior in the end. For instance, it is not clear 
beforehand, whether the approximation of infinite size (length, area) of the 
linear or planar detectors works better than the zero dimension approxima- 
tion for point detectors. Further developments will resolve these questions. 

8.1.2 Direct imaging techniques 

Some direct imaging techniques have been suggested, which might not require 
mathematical reconstructions. See, for instance, [79J about an acoustic lens 
system. 

8.1.3 Using contrast agents 

Contrast agents to improve TAT imaging have been developed (e.g., [24J). 
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8.1.4 Passive thermoacoustic imaging 

The TAT model we have considered can be called "active thermoacoustic to- 
mography," due to the set-up when the practitioner creates the signal. There 
has been some recent development of the "passive thermoacoustic tomogra- 
phy," where the thermoacoustic signal is used to image the temperature 
sources present inside the body. One can find a survey of this area in [94] . 

8.2 Uniqueness 

8.2.1 Sketch of the proof of Theorem [TO] 

We provide here a brief outline of the rather technical proof of Theorem [TUJ 
Suppose that / is compactly supported, not identically zero, and such 
that Rgf = 0. Our previous considerations show that one can assume that 
S is an algebraic curve (not a straight line) that is contained in the set of 
zeros of a non-trivial harmonic polynomial. Now one touches the boundary 
of the support of / from outside by a circle centered on S. Then microlocal 
analysis of the operator R$ (which happens to be an analytic Fourier Integral 
Operator, FIO pH 021 SHI S3 113 ESI [98]) shows that, due to the equality 
Rsf = 0, at the tangency point the vector co-normal to the sphere should not 
belong to the analytic wave front of / (microlocal regularity of solutions of 
Rsf = 0). This, for instance, can be also extracted from the results of [105] . 
On the other hand, a theorem by Hormander and Kashiwara [5H Theorem 
8.5.6] shows that this vector must be in the analytic wave front set, since 
/ = on one side of the sphere (a microlocal version of uniqueness of analytic 
continuation). This way, one gets a contradiction. Unfortunately, the life is 
not so easy, and the proof sketched above does not go through smoothly, due 
to possible cancelation of wavefronts at different tangency points. Then one 
has to involve the geometry of zeros of harmonic polynomials [37] to exclude 
the possibility of such a cancelation. 

Thus, the proof uses microlocal analysis and geometry of zeros of har- 
monic polynomials. Both these tools have their limitations. For instance, 
the microlocal approach (at least, in the form it is used in [7J) does not al- 
low considerations of non-compactly supported functions. Thus, the validity 
of the Theorem for arbitrarily fast decaying, but not compactly supported, 
functions is still not established, albeit it most certainly holds. On the other 
hand, the geometric part does not work that well in dimensions higher than 
two. Development of new approaches is apparently needed in order to over- 
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come these hurdles. A much simpler PDE approach has emerged recently 
[33] (see also [T2] and the next Section), albeit its achievements have been 
limited so far. 

8.2.2 Some open problems concerning uniqueness 

As it has already been mentioned, one can consider the practical problems 
about uniqueness resolved. However, the mathematical understanding of the 
uniqueness problem for the restricted spherical mean operators R$ is still 
unsatisfactory. Here are some questions that still await their resolution: 

1. Describe uniqueness sets in dimensions larger than 2 (prove the Con- 
jecture SJ). Recent limited progress, as well as variations on this theme 
can be found in [T]-[T2]. 

2. Prove Theorem [3] without using microlocal and harmonic polynomial 
tools. 

3. Prove Theorem [3] on uniqueness sets S under the condition of suffi- 
ciently fast decay (rather than compactness of support) of the function. 
Very little is known for the case of functions without compact support. 
The main known result is of [3], which describes for which values of 
1 < p < oo the result of Corollary [2] still holds: 

Theorem 12 f3j/ Let S be the boundary of a bounded domain in M. d and 
f e L p (R d ) such that R s f = 0. Ifp < 2d/(d-l), then f = (and thus 
S is injectivity set for this space). This fails for any p > 2d/ (d — 1). 

8.3 Inversion 

Albeit closed form (backprojection type) inversion formulas are available now 
for the cases of S being a plane (and object on one side from it), cylinder, 
and a sphere, there is still some mystery surrounding this issue. 

1. Can one write a backprojection type inversion formula in the case of 
the constant sound speed for a closed surface S which is not a sphere? 
We suspect that the answer to this question is negative (see also related 
discussion in [T5| [27]). 
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2. The inversion formulas for S being a sphere assume that the object to 
be imaged is inside S. One can check on simplest examples that if the 
support of function f(x) reaches outside S, the inversion formulas do 
not reconstruct the function correctly even inside of S. See [5] for a 
discussion. 

3. The I. Gelfand's school of integral geometry has developed a marvelous 
machinery of the so called k operator, which provides a general ap- 
proach to inversion and range descriptions for transforms of Radon type 
[3"%1 13"§] . In particular, it has been applied to the case of integration of 
various collections ( "complexes" ) of spheres in [3"5| |4"T] . This consider- 
ation seems to suggest that one should not expect explicit closed form 
inversion formulas for R$ when S is a sphere. We, however, know that 
such formulas have been discovered recently [33j 69J. This apparent 
controversy has not been resolved. 

4. Can one derive any more explicit analytic formulas from (128]) ? 

5. Can the series expansion formulas of Theorem [8] be efficiently imple- 
mented? 

One can also mention that in some works [151 123] it is suggested to use in 
the TAT problem not only the values of the pressure measured by transducers 
on the observation surface S, but its normal derivative to S as well. If one 
knows both, then taking Fourier transform in the time variable and using the 
whole space Green's function for the Helmholtz equation leads immediately 
to a reconstruction formula for the solution (which seems to be much simpler 
than what is proposed in [23] ). The problem is that this normal derivative is 
not measured by TAT devices. Under some circumstances (e.g., when there 
are no sources of ultrasound outside S), one can prove the theoretical pos- 
sibility of recovering the missing normal derivative. This, however, does not 
seem to us to be a plausible procedure. In rare cases (planar, cylindrical, 
or spherical surface S), when involvement of the normal derivative can be 
eliminated (e.g., [27]), this might lead to feasible inversion algorithms, 
but in these explained before in this text, explicit and nicely im- 

plementable analytic inversion formulas are available. So, jury is still out on 
this issue as well. 
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8.4 Stability 

Stability of inversion when S is a sphere surrounding the support of f(x) is 
the same as for the standard Radon transform, as the results of [91] and sec- 
ond statement of Theorem [TT1 show. However, if the support reaches outside, 
albeit Corollary [2] still guarantees uniqueness of reconstruction, stability (at 
least for the parts outside S) is gone. Indeed, Theorem [9] shows that some 
parts of singularities of / outside S will not be stably "visible." 

8.5 Range 

As Theorem IH] states, the range conditions 2 and 3 of Theorem [TD] are neces- 
sary also for non-spherical closed surfaces S and for functions with support 
outside S. They, however, are not expected to be sufficient, since Theorem [9] 
indicates that one might expect non-closed ranges in some cases. The same 
applies for non-constant sound speed case. 
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